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Abstract This review discusses Smoothed Particle Hydrodynamics (SPH) in the astrophysical 
context, with a focus on inviscid gas dynamics. The particle-based SPH technique allows an 
intuitive and simple formulation of hydrodynamics that has excellent conservation properties 
and can be coupled to self-gravity easily and highly accurately. The Lagrangian character of 
SPH allows it to automatically adjust its resolution to the clumping of matter, a property that 
makes the scheme ideal for many applications in astrophysics, where often a large dynamic 
range in density is encountered. We discuss the derivation of the basic SPH equations in their 
modern formulation, and give an overview about extensions of SPH developed to treat physics 
such as radiative transfer, thermal conduction, relativistic dynamics or magnetic fields. We also 
briefly describe some of the most important applications areas of SPH in astrophysical research. 
Finally, we provide a critical discussion of the accuracy of SPH for different hydrodynamical 
problems, including measurements of its convergence rate for important classes of problems. 
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1 INTRODUCTION 

Smoothed Particle Hydrodynamics (SPH) is a technique for approximating the contin 
uum dynamics of fluids through the use of particles, which may also be viewed as inter- 



polatio n po ints. SPH was originally dev eloped in astrophysics, as introduced by Lucy 
( 19771) and iGingold & Monaghanl (1 19771) a good 30 years ago. Since then it has found 
widespread use also in other areas of science and engineering. In this review, I discuss 
SPH in its modern form, based on a formulation derived from variational principals, 
giving SPH very good conservative properties and making its derivation largely free of 

ad-hoc choices that needed to be made in older versions of SPH. 

A few excellent reviews have dis cussed SPH previously (e.g. Monagharl 19921 2005 



Dolag et all 2008; iRosswogl. 12009b . hence we will concentrate primarily on recent de 



velopments and on a critical discussion of SPH's advantages and disadvantages, rather 
than on giving a full historical account of the most important literature on SPH. Also, we 
will for the most part restrict the discussion of SPH to inviscid ideal gases, which is the 
relevant case for the most common applications in astrophysics, especially in cosmology. 
Only in passing we comment on some other important uses of SPH in astronomy where 
fluids quite different from an ideal gas are modelled, e.g. in planet formation. Fully 
outside the scope of this review are the many successful applications of SPH-based tech- 
niques in fields such as geophysics or engineering. For example, free-surface flows tend 
to be very difficult to model with Eulerian methods, while this is comparatively easy with 
SPH. As a result, there are many applications of SPH to problems such as dam-braking, 
avalanches, and the like, which we however will not discuss here. 

Numerical simulations have become an important tool in astrophysical research. For 
example, cosmological simulations of structure formation within the ACDM model 
model have been instrumental to understand the non-linear outcome of the initial condi- 
tions predicted by the theory of inflation. By now, techni ques to simulate collis i onless 



dark matter through the particle-based N-body method (IHockney & Eastwood 119811) 
have fully matured and are comparatively well understood. However, to represent the 
collisional baryons as well, the hydrodynamical fluid equations need to be solved, which 
represents a much harder problem than the dark matter dynamics. On top of the more 
complicated gasdynamics, additional physics, like radiation fields, magnetic fields, or 
non-thermal particle components need to be numerically followed as well to produce 
realistic models of the formation and evolution of galaxies, stars or planets. There is 
therefore ample need for robust, accurate, and efficient hydrodynamical discretization 
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techniques in astrophysics. 

The principal idea of SPH is to treat hydrodynamics in a completely mesh-free fash- 
ion, in terms of a set of sampling particles. Hydrodynamical equations of motion are then 
derived for these particles, yielding a quite simple and intuitive formulation of gas dy- 
namics. Moreover, it turns out that the particle representation of SPH has excellent con- 
servation properties. Energy, linear momentum, angular momentum, mass, and entropy 
(if no artificial viscosity operates) are all simultaneously conserved. In addition, there are 
no advection errors in SPH, and the scheme is fully Galilean invariant, unlike alternative 
mesh-based Eulerian techniques. Due to its Lagrangian character, the local resolution 
of SPH follows the mass flow automatically, a property that is extremely convenient 
in representing the large density contrasts often encountered in astrophysical problems. 
Together with the ease in which SPH can be combined with accurate treatments of self- 
gravity (suitable and efficient gravity solvers for particles can be conveniently taken from 
a cosmological N-body code developed for the representation of dark matter), this has 
made the method very popular for studying a wide array of problems in astrophysics, 
ranging from cosmological structure growth driven by gravitational instability to studies 
of the collisions of protoplanets. Furthermore, additional 'sub-resolution' treatments of 
unresolved physical processes (such as star formation in galaxy-scale simulations) can 
be intuitively added at the particle level in SPH. 

In this review, we first give, in Section|2j a derivation of what has become the standard 
formulation of SPH for ideal gases, including also a description of the role of artificial 
viscosity. We then summarize some extensions of SPH to include additional physics 
such as magnetic fields or thermal conduction in Section [3j followed by a brief descrip- 
tion of the different application areas of SPH in astronomy in Section H] Despite the 
popularity of SPH, there have been few systematic studies of the accuracy of SPH when 
compared with the traditional Eulerian approaches. We therefore include a discussion of 
the convergence, consistency and stability of standard SPH in Section [5j based on some 
of our own tests. Finally, we give a discussion of potential future directions of SPH 
development in Section [6l and our conclusions in Section [7] 



2 BASIC FORMULATION OF SPH FOR IDEAL GASES 



2.1 Kernel Interpolants 

At the heart of smoothed p article hydrodynamics l ie so-c alled kernel interpolants, which 
are discussed in detail by iGingold & Monaghanl dl982l) . In particular, we use a kernel 
summation interpolant for estimating the density, which then determines the rest of the 
basic SPH equations through the variational formalism. 

For any field F(r), we may define a smoothed interpolated version, F s (r), through a 
convolution with a kernel W(r,h): 

F,(r) = J F(r)W(r-r',h)dr'. (1) 

Here h describes the characteristic width of the kernel, which is normalized to unity 
and approximates a Dirac 8-function in the limit h — > 0. We further require that the 
kernel is symmetric and sufficiently smooth to make it at least differentiable twice. On e 



possibility for W is a Gaussian, which was in fact used by IGingold & Monaghanl (I l977|) 



However, most current SPH implementations are based on kernels with a finite support. 
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Usually a cubic spline is adopted with W(r,h) = w(£), and 

' \-6q 2 + 6q\ 0<q< L 

W3d(<?) = - < 



71 



2(1 -qf, \<q<l, (2) 
0, q>l, 



in three-dimensional normal ization. This kernel belongs to a broader class of interpola- 
tion a nd smoothing kernels ( Schoenberglll969l : lHockney & Eastwoodlll98~il lMonaghan. 



19851) . Note that in the above most commonly used definition of the smoothing length 
h, the kernel drops to zero at a distance of r = 2h. Through Taylor expansion, it is easy 
to see that the kernel interpolant is at least second-order accurate due to the symmetry of 
the kernel. 

Suppose now we know the field at a set of points r ; , i.e. Fj = Ffa). The points have 
an associated mass m ; and density p,-, such that Ar, ~ m,-/p,- is their associated finite 
volume element. Provided the points sufficiently densely sample the kernel volume, we 
can approximate the integral in Eqn. (OQ) with the sum 

Fs(r)-L^ F jW(r-rj,h). (3) 
; PJ 

This is effectively a Monte-Carlo integration, except that thanks to the comparatively 
regular distribution of points encountered in practice, the accuracy is considerably better 
than for a random distribution of the sampling points. In particular, for points in one di- 
mension with equal spacing d, one can show that for h = d the sum of Eqn. (O provides 
a second order accurate approximation to the real underlying function. Unfortunately, 
for the irregular yet somewhat ordered particle configurations encountered in real appli- 
cations, a formal error analysis is not straightforward. It is clear however, that at the very 
least one should have h>d, which translates to a minimum of ~ 33 neighbours in 3D. 

Importantly, we see that the estimate for F s (r) is defined everywhere (not only at 
the underlying points), and is differentiable thanks to the differentiability of the kernel, 
albeit with a considerably higher interpolation error for the derivative. Moreover, if we 
set F(r) = p(r), we obtain 

p s (r)~Y f mjW(r-r j ,h) 1 (4) 

j 

yielding a density estimate just based on the particle coordinates and their masses. In 
general, the smoothing length can be made variable in space, h = h(r,t), to account 
for variations in the sampling density. This adaptivity is one of the key advantages 
of SPH and is essentially always used in practice. There are two options to introduce 
the variability of h into Eqn. ((U). One is by adopting W(r — Tj, h(r)) as kernel, which 
corresponds to the 'scatter' approach (IHernquist & K atz. 1989). It has the advantage 



that the volume integral of the smoothed field recovers the total mass, / p s (r) dr = m,-. 
On the other hand, the so-called 'gather' approach, where we use W(r — ry,/i(r,-)) as 
kernel in Eqn. (01), requires only knowledge of the smoothing length hi for estimating 
the density of particle i, which leads to computationally convenient expressions when 
the variation of the smoothing length is consistently included in the SPH equations of 
motion. Since the density is only needed at the coordinates of the particles and the total 
mass is conserved anyway (since it is tied to the particles), it is not important that the 
volume integral of the gather form of p. v (r) exactly equals the total mass. 

In the following we drop the subscript s for indicating the smoothed field, and adopt 
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as SPH estimate of the density of particle i the expression 

N 

p i = Y t m J W(T i -Tj,hi). (5) 

7=1 

It is clear now why kernels with a finite support are preferred. They allow the summation 
to be restricted to the Af ngb neighbors that lie within the spherical region of radius 2h 
around the target point r,-, corresponding to a computational cost of order o(N ng bN) for 
the full density estimate. Normally this number N ng b of neighbors within the support of 
the kernel is approximately (or exactly) kept constant by choosing the hi appropriately. 
/V ng b hence represents an important parameter of the SPH method and needs to be made 
large enough to provide sufficient sampling of the kernel volumes. Kernels like the 
Gaussian on the other hand would require a summation over all particles N for every 
target particle, resulting in a O (N 2 ) scaling of the computational cost. 

If SPH was really a Monte-Carlo method, the accuracy expected from the interpola- 
tion errors of the density estimate would be rather problematic. But the errors are much 
smaller because the particles do not sample the fluid in a Poissonian fashion. Instead, 
their distances tend to equilibrate due to the pressure forces, which makes the interpola- 
tion errors much smaller. Yet, they remain a significant source of error in SPH and are 
ultimately the primary origin of the noise inherent in SPH results. 

Even though we have based most of the above discussion on the density, the general 
kernel interpolation technique can also be applied to other fields, and to the construction 
of differential operators. For example, we may write down a smoothed velocity field and 
take its derivative to estimate the local velocity divergence, yielding: 

(V-v^E-^vy-V.-Wfc-ry,*). (6) 

j p j 

However, an alternative estimate can be obtained by considering the identity pV • v = 
V(pv) — v • Vp, and computing kernel estimates for the two terms on the right hand side 
independently. Their difference then yields 

(V • v)j = - £m;(vj - v«) • ViW(rj - r b h). (7) 

This pair-wise formulation turns out to be more accurate in practice. In particular, it 
has the advantage of always providing a vanishing velocity divergence if all particle 
velocities are equal. 

2.2 Variational Derivation 

The Euler equations for inviscid gas dynamics in Lagrangian (comoving) form are given 
by 

^7 + pV-v = 0, (8) 
at 

dv VP „ 
at p 
du P 

+ W = 0, (10) 
at p 
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where d/dt = d/dt + v • V is the convective derivative. This system o f partial differential 
equations expresses conservation of mass, momentum and energy. lEckartl (i960) has 



shown that the Euler equations for an inviscid ideal gas follow from the Lagrangian 

p(y-uW (11) 

This opens up an interesting route for obtaining discretized equations of motion for 
gas dynamics. Instead of working with the continuum equations directly and trying to 
heuristically work out a set of accurate difference formulas, one can discretize the La- 
grangian and then derive SPH equations of motion by applying the variational prin cipals 



of classical mechanics, an approach first proposed by Gingold & Monaghanl ( 1982|) . Us- 
ing a Lagrangian also immediately guarantees certain conservation laws and retains the 
geometric structure imposed by Hamiltonian dynamics on phase space. 

We here follow this el egant idea, whic h was first worked out bv Springel & Hern- 
quist ( 2002|) with a consistent accounting of variable smoothing lengths. We start by 



discretizing the Lagrangian in terms of fluid particles of mass m,-, yielding 

^SPH = £ Q m * V / ~ m ' U ^j > ( 12 ) 

where it has been assumed that the thermal energy per unit mass of a particle can be 
expressed through an entropic function A, of the particle, which simply labels its specific 
thermodynamic entropy. The pressure of the particles is 

Pi=A,pl=(y-l)piUi, (13) 

where yis the adiabatic index. Note that for isentropic flow (i.e. in the absence of shocks, 
and without mixing or thermal conduction) we expect the A, to be constant. We hence 
define m ( , the thermal energy per unit mass, in terms of the density estimate as 

u i (pt)=A i *- r[ . (14) 

This raises the question of how the smoothing lengths /z, needed for estimating p, should 
be determined. As we discussed above, we would like to ensure adaptive kernel sizes, 
meaning that the number of points in the kernel should be approximately constant. In 
much of the older SPH literature, the number of neighbours was allowed to vary within 
some (small) range around a target number. Sometimes the smoothing length itself was 
evolved with a differential equation in time, exploiting th e continuity relation and th e 



expectation that ph 3 should be approximately constant (e.g lSteinmetz & MuellerLll993h . 
In case the number of neighbours outside t he kernel happened to fall ou tside the allowed 
range, h was suitably readjusted. However. iNelson & Papaloizoul d 19941) pointed out that 



for smoothing lengths varied in this way, the energy is not conserved correctly. They 
showed that the errors could be made smaller by keeping the number of neighbours 
exactly constant, and they also derived leading order correction terms (which became 
known as Vh terms) for the classic SPH equations of motion that could reduce them still 
further. In the modern formulation discussed below, these V/z-terms do not occur; they 
are implicitly included at all orders. 

The central trick making this possible is to require that the mass in the kernel volume 
should be constant, viz. 

Pih] = const (15) 
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for three dimensions. Since p,- = p,-(ri,r2, . . .r#,/i,-) is only a function of the particle 
coordinates and of hj, this equation implicitly defines the function h t = h i (r l ,r 2 , ■ ■ .r N ) 
in terms of the particle coordinates. 

We can then proceed to derive the equations of motion from 



This first gives 



d dL dL _ 
At 3r,- 3r, 



dVi " Pj d Pj 

At ; fi Pj dr i 



(16) 



(17) 



where the derivative dpy/dr, stands for the total variation of the density with respect to 
the coordinate r,, including any variation of hj this may entail. We can hence write 



(18) 



where the smoothing length is kept constant in the first derivative on the right hand 
side (in our notation, the Nabla operator V; = d/drj means differentiation with respect 
to r, holding the smoothing lengths constant). On the other hand, differentiation of 
Pjhj = const with respect to r; yields 



dhj 3r, 



1 + 



3 Pj fdpj 



dhj 



-l 



(19) 



Combining equations (fT8T ) and ( fT9l > we then find 



dpj 



hj dpj 



dr, V 3p 7 dhj 



-l 



Using 



k=\ 



we finally obtain the equations of motion 



QI 7=1 



fAViWijM + fAViWijihj) 
Pi Pj 



where the /; are defined by 



hi_dpi 
h 3p, 01, 



(20) 



(21) 



(22) 



(23) 



and the abbreviation Wjj{h) = W(\ri — rj\,h) has been used. Note that the correction 
factors can be easily calculated alongside the density estimate, all that is required 
is an additional summation to get 3p,/3r,- for each particle. This quantity is in fact 
also useful to get the correct smoothing radii by iteratively solving p,-/^ = const with a 
Newton-Raphson iteration. 

The equations of motion (1221 for inviscid hydrodynamics are remarkably simple. In 
essence, we have transformed a complicated system of partial differential equations into 
a much simpler set of ordinary differential equations. Furthermore, we only have to 



8 



V. Springel 



solve the momentum equation explicitly. The mass conservation equation as well as the 
total energy equation (and hence the thermal energy equation) are already taken care of, 
because the particle masses and their specific entropies stay constant for reversible gas 
dynamics. However, later we will introduce an artificial viscosity that is needed to allow 
a treatment of shocks. This will introduce additional terms in the equation of motion 
and requires the time integration of one thermodynamic quantity per particl e, which can 



qua ntity per partu 

either be chosen as entropy or thermal energy. Indeed, iMonaghanl (|2002|) pointed out 



that the above formulation can also be equivalently expressed in terms of thermal energy 
instead of entropy. This follows by taking the time derivative of Eqn. (fT4l ). which first 
yields 

^ = *J> r jj* (24) 



Using equations (1201) and (1211) then gives the evolution of the thermal energy as 

7*7 =fi^!Lrn ] {y i -y j )-VW ij {h i ), (25) 

which needs to be integrated along the equation of motion if one wants to use the ther- 
mal energy as independent thermodynamic variable. There is no difference however 
to using the entropy; the two are completely equi valent in the va riational formulation. 



This also solves the old problem pointed out by iHernquistl (119931) that the classic SPH 
equations did not properly conserve energy when the entropy was integrated, and vice 
versa. Arguably, it is numerically advantageous to integrate the entropy though, as this 
is computationally cheaper and eliminates time integration errors in solving Eqn. d25l ). 

Note that the above formulation readily fulfills the conservation laws of energy, mo- 
mentum and angular momentum. This can be shown based on the discretized form of the 
equations, but it is also manifest due to the symmetries of the Lagrangian that was used 
as a starting point. The absence of an explicit time dependence gives the energy conser- 
vation, the translational invariance implies momentum conservation, and the rotational 
invariance gives angular momentum conservation. 

Other derivations of the SPH equations based on constructing kernel interpolated ver- 
sions of differential operators and app lying them directly to the Euler equations are also 



possible (see, e.g.. IMonaghanl 11992). However, these derivations of 'classic SPH' are 
not unique in the sense that one is left with several different possibilities for the equations 
and certain ad-hoc symmetrizations need to be introduced. The choice for a particular 
formulation then n eeds to rely on exp erimentally comparing the performance of many 



different variants (Thacke r et all 12000). 



2.3 Artificial Viscosity 

Even when starting from perfectly smooth initial conditions, the gas dynamics described 
by the Euler equations may readil y produce true discontinu ities in the form of shock 



waves and contact discontinuities ( Landau & Lifshitzl Il959h . At such fronts the differ- 



ential form of the Euler equations breaks down, and their integral form (equivalent to the 
conservation laws) needs to be used. At a shock front, this yields the Rankine-Hugoniot 
jump conditions that relate the upstream and downstream states of the fluid. These rela- 
tions show that the specific entropy of the gas always increases at a shock front, implying 
that in the shock layer itself the gas dynamics can no longer be described as inviscid. In 
turn, this also implies that the discretized SPH equations we derived above can not cor- 
rectly describe a shock, simply because they keep the entropy strictly constant. 
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One thus must allow for a modification of the dynamics at shocks and somehow in- 
troduce the necessary dissipation. This is usually accomplished in SPH by an artificial 
viscosity. Its purpose is to dissipate kinetic energy into heat and to produce entropy in 
the process. The usual approach is to parameterize the artificial viscosity in terms of 
a friction force that damps the relative motion of particles. Through the viscosity, the 
shock is broadened into a resolvable layer, something that makes a description of the 
dynamics everywhere in terms of the differential form possible. It may seem a daunting 
task though to somehow tune the strength of the artificial viscosity such that just the 
right amount of entropy is generated in a shock. Fortunately, this is however relatively 
unproblematic. Provided the viscosity is introduced into the dynamics in a conservative 
fashion, the conservation laws themselves ensure that the right amount of dissipation 
occurs at a shock front. 

What is more problematic is to devise the viscosity such that it is only active when 
there is really a shock present. If it also operates outside of shocks, even if only at a 
weak level, the dynamics may begin to deviate from that of an ideal gas. 

The viscous force is most often added to the equation of motion as 



dv. 
At 



= -J>,-nyViWy, 

vise j= 1 



(26) 



where 



_ 1 



^[Wijihij+Wijihj)] 



(27) 



'til- 



denotes a symmetrized kernel, which some researchers prefer to define as Wjj = Wjj 
hj]/2). Provided the viscosity factor Ily is symmetric in i and j, the viscous force be- 
tween any pair of interacting particles will be antisymmetric and along the line joining 
the particles. Hence linear momentum and angular momentum are still preserved. In 
order to conserve total energy, we need to compensate the work done against the viscous 
force in the thermal reservoir, described either in terms of entropy, 



dAi 
~d7 



1 y-1 



N 



2 Pi i=l 



(28) 



or in terms of thermal energy per unit mass, 



~dt 



X>,n, 



(29) 



'7=1 



where v,-y = v,- — Vj. There is substantial freedom in the detailed parametrization of 
the viscosity IT, ; . The most commonly used formula tion is an improved version of the 
viscosity introduced by lMonaghan & Gingoldl (119830 . 



n, 



-CLCijWj + (3//r. /Pij 



ifv l7 -r l7 <0 
otherwise, 



(30) 



with 



(31) 



Here hjj and p,y denote arithmetic means of the corresponding quantities for the two 
particles i and j, with Cjj giving the mean sound speed, whereas r,y = r, — ry. The 



10 



V. Springel 



strength of the viscosity is regulated by the parameters a and p\ with typical values in 
the range a ~ 0.5 — 1.0 and the frequent choice of [3 = 2a. The parameter £ ~ 0.01 is 
introduced to protect against singularities if two particles happen to get very close. 

In this form, the artificial viscosity is basically a combination of a bulk and a von 
Neumann-Richtmyer viscosity. Historically, the quadratic term in fiy has been added to 
the original Monaghan-Gingold form to prevent particle penetration in high Mach num- 
ber shocks. Note that the viscosity only acts for particles that rapidly approach each 
other, hence the entropy production is always positive definite. Also, the viscosity van- 
ishes f or solid- b ody ro tation, but not for pure shear flows. To cure this problem in shear 



flows, iBalsaral (1 1995b suggested adding a correction factor to the viscosity, reducing 
its strength when the shear is strong. This can be achieved by multiplying Tljj with a 
prefactor (ff^ + /j w )/2, where the factors 

f> AV = iv T/lv l < 32 > 

|V-v|; + |V X v|; 

are meant to measure the rate of local compression in relation to the strength of the local 
shear (estimated with formulas such as Eqn. |7]). This Balsara switch has often been 
successfully used in multi-dimensional flows and is enabled as default in many SPH 
codes. We note however that it may be problematic sometimes in cases where shocks 
and shear occur together, e.g. in oblique shocks in differentially rotating disks. 

In some studies, al ternativ e forms of viscosity have b een test ed. For example. Mon- 
aghan d 19971) proposed a modified form of the viscosity based on an analogy to the 
Riemann problem, which can be written as 

2 pij 

where v-j g = [c, + cj — 3w;y] is an estimate of the signal velocity between two particles i 
and j, and = Vy • ry/ ry | is the relative velocity projected onto the separation vector. 
This viscosity is identical to (l30l if one sets [3 = 3/2 a and replaces Wjj with wj. The 
main difference between the two viscosities lies therefore in the additional factor of 
hij/rij that mj carries with respect to Wij. In equations (l30l and (|3"TT ). this factor weights 
the viscous force towards particle pairs with small separations. In fact, after multiplying 
with the kernel derivative, this weighting is strong enough to make the viscous force in 
equation ( TSTT ) diverge as <=c 1 / ry for small pair separations up to the limit set by the s/j? 
term. 



Lombardi et al. (1999) have systematically tested different parametrization of viscos- 



ity, but in general, the standard form was found to w ork best. Recently the use of a tensor 
artificial viscosity was conjectured by lQwenl ( 2004 ) as part of an attempt to optimize the 



spatial resolution of SPH. However, a disadvantage of this scheme is that it breaks the 
strict conservation of angular momentum. 

In attempting to reduce the numerical viscosity of SPH in regions away from shocks, 
several studies have instead advanced the idea of keeping the functional form of the 
artificial viscosity, but making the viscosity s trength parameter a variable in time. Such a 



scheme was first suggested by lMorrisl ( 19971) . and was successfully appl ied for studying 



astrophysical turbule nce more faithfu lly in SPH dDolag et all [2005b), and to follow 



neutron star mergers (IRossw og. 2005). Adopting [3 = 2a, one may evolve the parameter 



a individually for each particle with an equation such as 
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where Sj is some source function meant to ramp up the viscosity rapidly if a shock 
is detected, while the first term lets the viscosity exponentially decay again to a pre- 
scribed minimum value a m j n on a timescale x,. So far, simple source functions like 
Si = max[— (V • v),-,0] and timescales x, ~ hj/ct have been explored and the viscosity a, 
has often also been prevented from becoming higher than some prescribed maximum 
value a max . It is clear that the success of such a variable a scheme depends critically 
on an appropriate source function. The form above can still not distinguish purely adi- 
abatic compression from that in a shock, so is not completely free of creating unwanted 
viscosity. 



2.4 Coupling to Self-gravity 

Self-gravity is extremely important in many astrophysical flows, quite in contrast to other 
areas of computational fluid dynamics, where only external gravitational fields play a 
role. It is noteworthy that Eulerian me sh-based approaches do not manifest l y co nserve 
total energy if self-gravity is included ( Miiller & Steinmetz , 19951: ISpringelL 2009|) . but 



it can be easily and accurately incorporated in SPH. 

Arguably the best approach to account for gravity is to include it directly into the dis- 
cretized SPH Lagrangian, which has the advantage of also allowing a consistent treat- 
ment of variable gravitational softening lengths. This was first con ducted by Price & 
Monaghan d2007l) i n the context of adaptiv e resolution N-body methods for collisionless 



dynamics (see also lBagla & Khandai , 120091) . 



Let <J>(r) = G^, m ! (j)(r — r ,-,£,-) be the gravitational field described by the SPH point 
set, where £, is the gravitational softening length of particle i. We then define the total 
gravitational self-energy of the system of SPH particles as 

1 G 
£ P ot = ~Y, m M r i) = T^m/ffl^^.Ej). (35) 
L i L ij 

For SPH including self-gravity, the Lagrangian then becomes 

LspH = 11 (^mrf ~ mu^j - ^Y^mimj^(rij,Zj). (36) 

As a result, the equation of motion acquires an additional term due to the gravitational 
forces, given by: 



nii a 



;rav dE pot 



dri 
Y^Gniiinj 



Tij W(rij,Et)+tf(rtj,Ej)] 



J 



2 



-2^°™^^^^ (37) 

where (j)'(r,£) = 3(j)/3r. 

The first sum on the right hand side describes the ordinary gravitational force, where 
the interaction is symmetrized by averaging the forces in case the softening lengths be- 
tween an interacting pair are different. The second sum gives an additional force com- 
ponent that arises when the gravitational softening lengths are a function of the particle 
coordinates themselves, i.e. if one invokes adaptive gravitational softening. This term 
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has to be included to make the system properly conservative wh en spatially adaptive 
gravitational softening lengths are used (|Price & MonaghanL 120071) . 

In most cosmological SPH codes of structure formation, this has usually not been 
done thus far and a fixed gravitational softening was used for collisionless dark matter 
particles and SPH particles alike. However, especially in the context of gravitational 
fragmentation problems , it appears natural and also indicated by accuracy considerations 
( Bate & Burkerti 1997 ) to tie the gravitational softening length to the SPH smoothing 
length, even though some studies also caution against this strategy (Williams, Churches 
& Nelson. 120041) . Adopting £, = hi, and determining the smoothing lengths as described 
above, we can calculate the 3e//3r; = dhj/drj term by means of equations ( fl9l ) and (f2T]). 
Defining the quantities 



hj ^ 3<|>(r 



1j 



where the factors /, are defined as in Eqn. ((23 
acceleration compactly as 



jk,hj) 



dh 



(38) 



we can then write the gravitational 



dv; 
At 



Yij W(rij,ht)+tf(rij,hj) 



j r 'J 



G 



+ ~L m J fo/V/WyCftj) +i\jViWij{hj)] . 



(39) 



This acceleration has to be added to the ordinary SPH equation of motion (1221 ) due to 
the pressure forces, which arises from the first part of the Lagrangian (l36l) . Note that 
the factors r\j can be conveniently calculated alongside the SPH density estimate. The 
calculation of the gravitational correction force, i.e. the second sum in Eqn. (|39l >, can 
then be done together with the usual calculation of the SPH pressure forces. Unlike in 
the primary gravitational force, here only the nearest neighbours contribute. 

We note that it appears natural to relate the functional form of the gravitational soften- 
ing to the shape of the SPH smoothing kernel, even though this is not strictly necessary. 
If this is done, §(r,h) is determined as solution of Poisson's equation, 



V 2 ty(r,h) =4nW(r,h). 



(40) 



The explicit functional form for §(r,h) resulting from this identification for the cubic 
spline kernel is also often e mployed in collisionless N-bqdy cod es and can be found, for 
example, in the appendix of iSpringel. Yoshida & White! (|2001a|) . 



2.5 Implementation Aspects 



The actual use of the discretized SPH equations in simulation models requires a time- 
integration scheme. The Hamil tonian character of SPH allows in principle the u se of 
symplectic integration schemes (Hai rer. Lubich & Wanneii 1200 2; Springe 1, 2005) such 
as the leapfrog, which respects the geometric phase-space constraints imposed by the 
conservation laws and can prevent the build-up of secular integration errors with time. 
However, since most hydrodynamical problems of interest are not reversible anyway, 
this aspect is less important than in dissipation-free collisionless dynamics. Hence any 
second-order accurate time integration scheme, like a simple Runge-Kutta or predictor- 
corrector scheme, may equally well be used. It is common practice to use individual 
timesteps for the SPH particles, often in a bloc k-structured scheme wi th particles ar- 
ranged in a power-of-two hierarchy of timesteps (Her nquist & Katzi Il989h . This greatly 
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increases the efficiency of calculations in systems with large dynamic range in timescale, 
but can make the optimum choice of timesteps tricky. For example, the standard Courant 
timestep criterion of the form 

A?;=C CFL -, (41) 

Ci 

usually used in SPH, where Ccfl ~ 0.1 —0.3 is a dimensionless parameter, may not 
guarantee fine enough time-stepping ahead of a blast wave propagating into very cold 
gas. This problem ca n be avoided with an imp roved method for determining the sizes of 



individual timesteps dSaitoh & Makinol 12009). 



If self-gravity is included, one can draw from the algorithms employed in collisionless 
N-body dynamics to calculate the gravitational forces efficiently, such as tree-methods 
(Barnes or mesh-based gravity solvers. The tree approach, which provides 



for a hierarchical grouping of the particles, can also be used to address the primary 
algorithmic requirement to write an efficient SPH code, namely the need to efficiently 
find the neighboring particles inside the smoothing kernel of a particle. If this is done 
naively, by computing the distance to all other particles, a prohibitive o(N 2 ) scaling 
of the computational cost results. Using a range-searching technique together with the 
tree, this cost can be reduced to o(N ng \,N\ogN), independent of the particle clustering. 
To this end, a special walk of the tree is carried out for neighbor searching in which 
a tree node is only opened if there is a spatial overlap between the tree node and the 
smoothing kernel of the target particle, otherwise the corresponding branch of the tree 
can be immediately discarded. 

These convenient properties of the tree algorithm have been exploited for the develop- 
ment of a number of efficient SPH codes in astronomy, several of them are publicly avail- 
able and parallelized for distributed memory machines. Amons them are TreeSPH (Hern- 
quist & Katz. ll989l:lKatz. Weinberg & Hernquistl.ll996h. H YD RA dPearce & Couchman , 



et al.. l2009h . 



1997b . GADGET dSprineel. Y oshida & White, 2001b; Springel, 2005), GASOLINE (Wad- 
sley, Stadel & Ouinn. l2004h. MAGMA dRosswog & PriceL l2007h. and VINE (" Wetzstein 



3 EXTENSIONS OF SPH 



For many astrophysical applications, additional physical processes in the gas phase be- 
sides inviscid hydrodynamics need to be modelled. This includes, for example, magnetic 
fields, transport processes such as thermal conduction or physical viscosity, or radiative 
transfer. SPH also needs to be modified for the treatment of fluids that move at relativis- 
tic speeds or in relativistically deep potentials. Below we give a brief overview of some 
of the extensions of SPH that have been developed to study this physics. 



3.1 Magnetic Fields 

Magnetic fields are ubiquitous in astrophysical plasmas, where the conductivity can of- 
ten be approximated as being effectively infinite. In this limit one aims to simulate ideal, 
non-resistive magneto-hydrodynamics (MHD), which is thought to be potentially very 
important in many situations, in particular in star formation, cosmological structure for- 
mation, or accretion disks. The equations of ideal MHD are comprised of the induction 
equation 

^ = (B-V)v-B(V-v) (42) 



14 



V. Springel 



and the magnetic Lorentz force. The latter can be obtained from the Maxwell stress 
tensor 

1 / „ 1 



9vl 



'J 



as 



4* » BfB ' 



d<M; 



(43) 



(44) 



Working with the stress tensor is advantageous for deriving equations of motion that are 
discretized in a symmetric fashion. The magnetic force F then has to be added to the 
usual forces from gas pressure and the gravitational field. 

Firs t implementations of magnetic forces in SPH date back to Gingold & Monaghan 



d 19771) . soon followed by full implementations of MHD in SPH (Phillips & Monaghan, 
1985h . However, a significant problem with MHD in SPH has become apparent early 
on. The constraint equation V • B = 0, which is maintained by the continuum form 
of the ideal MHD equations, is in general not preserved by discretized versions of the 
equations. Those tend to build up some V • B ^ error over time, corresponding to 
unphysical magnetic monopole sources. In contrast, in the most modern mesh-based 
approaches to MHD, so-called constrained transport schemes (|Evans & Hawleyl Il988|) 
are able to accurately evolve the discretized magnetic field while keeping a vanishing 
divergence of the field. 

Much of the recent research on developing improved SPH realizations of MHD has 
therefore concentrated on constructing form ulations that either eliminate t he V • B^O 
error, or at l east keep it reasonably small dDolag. Ba rtelman n & Leschl 19991: Price 



& Monaghan, 2004; 



eep 
lab, 



120051 : IPricei 120091: iDolag & Stasvszvnl l2009t) . To this extent, 
different approaches hav e been investigated in the literature, involving periodic field 
cleaning techniques (e.g. lDedner et aUl2002l) . formulations of the equations that 'diffuse 
away' the terms if they should arise, or the use of the Euler potentials or the 

vector potential. The latter may seem like the most obvious solution, since deriving 
the magnetic field as B = V x A from the ve ctor potential A will automatically create a 
divergence free field. However, IPricd d2009h explored the use of the vector potential in 
SPH and concluded that there are substantial instabilities in this approach, rendering it 
essentially useless in practice. 

Another seemingly clean solution lies in the so called-Euler potentials (also known as 
Clebsch variables). One may write the magnetic field as the cross product B = Va x Vf3 
of the gradients of two scalar fields a and p\ In ideal MHD, where the magnetic flux is 
locked into the flow, one obtains the correct field evolution by simply advecting these 
Euler potentials a and f3 with the motion of the gas. This suggests a deceptively simple 
approach to ideal MHD; just construct the fields a and [3 for a given magnetic field, 
and then move these scalars along with the gas. Although this has been shown to work 
reasonably well in a number of simple test problems and also ha s been used in some real- 



world applications (Pri ce & BateL 120071 : iKotarba et all 120091) . it likely has significant 



problems in general MHD dynamics, as a result of the noise in SPH an d the inability 
of this scheme to account for any magnetic reconnection. Furthermore, Brandenburg! 

ooints out that the use of Euler potentials does not converge to a proper solution 
in hydromagnetic turbulence. 

Present formulations of MHD in SPH are hence back to d iscretizing the classic mag- 
netic induction equation. For example, IDolag & Stasvszvnl ^2009) give an implemen- 
tation of SPH in the GADGET code where the magnetic induction equation is adopted 
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as 



At p; 



J>; {B,[v y • ViWijfc)] - v i7 [Bi • ViWijihi)]} 



and the acceleration due to magnetic forces as 



— — = > mi 
At 4- y 



(45) 



(46) 



where the /)• are the correction factors of Eqn. d23]) , and Mj is the stress tensor of particle 
i. Combined with field cleaning techniques (|BOrve. Omang & TrulsenL 120011) and artifi- 
cial magnetic dissipation to keep V • B ^ errors under control, this leads to quite accu- 
rate results for many standard tests of MHD, like magnetic shock tubes or the Orszang- 
Tang vortex. A similar implementation of MHD in SPH is given in the independent 
MAGMA code bv lRosswog & Pried d2007l) . 



3.2 Thermal Conduction 



Diffusion processes governed by variants of the equation 



d<2 o 



(47) 



where Q is some conserved scalar field and D is a diffusion constant, require a discretiza- 
tion of the Laplace operator in SPH. This could be done in principle by differentiating 
a kernel interpolated version of Q twice. However, such a discretization turns out to be 
quite sensitive to the local particle distribution, or in other words, it is fai rly noisy. It 
is muc h better to use an approximation of the V 2 operator proposed first by iBrookshaw 
dl985h . in the form 

mi Qi — Qi „ — 

; ; - --V/Wy. 

J "' 



v 2 e| 



r 1 - 
'J 



(48) 



Clearv & Monaghanl d 19991) and ljubelgas. Springel & Dolad d2004 have used this to 
construct implementations of thermal conduction in SPH which allow for spatially vari- 
able conductivities. The heat conduction equation, 



where K is the heat conductivity, can then be discretized in SPH as follows, 



~d7 



in 



j [K, ■ Kj)(T; T, 



P»-Py 



Vii-ViW 



(49) 



(50) 



with Tjj = r ( - — Yj. In this form, the energy exchange between two particles is balanced 
on a pairwise basis, and heat always flows from higher to lower temperature. Also, it 
is easy to see that the total entropy increases in the process. This formulation has been 
used, for example, to study the influence of Spitzer conductivity d ue to electron transp ort 
on the thermal structure of the plasma in massive galaxy clusters dDolag et all 12004). 



3.3 Physical Viscosity 

If thermal conduction is not strongly suppressed by magnetic fields in hot plasmas, then 
one also expects a residual physical viscosity (a pure shear viscosity in this case). In 
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general, real gases may possess both physical bulk and physical shear viscosity. They 
are then correctly described by the Navier-Stokes equations, and not the Euler equations 
of inviscid gas dynamics that we have discussed thus far. The stress tensor of the gas can 
be written as 

' 3v; , dvj 2 c dvk \ y c dvk 



(51) 



where r\ and C, are the coefficients of shear and bulk viscosity, respectively. The Navier- 
Stokes equation including gravity is then given by 



dv VP „^ 1 

— = VO + -VG. 

dt p p 



(52) 



Sijacki & Springell (2006) have suggested an SPH discretization of this equation which 
estimates the shear viscosity of particle i based on 



dt 



ni; 



shear 



^ i [ v mj (h i )]\t + J ^[v i w ij (hj) 



pr 



(53) 



and the bulk viscosity as 



dv 
dt 



(,bulk j 



+ 



P} 



ViWijihj) 



(54) 



These viscous forces are antisymmetric, and together with a positive definite entropy 
evolution 



dM 
dt 



Y-l 



p; 



y-l 



2 pi ' p. 



(55) 



they conserve total energy. Here the derivatives of the velocity can be estimated with 
pair-wise formulations as in Eqn. f7]). 

In some sense, it is actually simpler for SPH to solve the Navier-Stokes equation than 
the Euler equations. As we discuss in more detail in Section 15.41 this is because the 
inclusion of an artificial viscosity in SPH makes the scheme problematic for purely in- 
viscid flow, since this typically introduces a certain level of spurious numerical viscosity 
also outside of shocks. If the gas has however a sizable amount of intrinsic physical 
viscosity anyway, it becomes much easier to correctly represent the fluid, provided the 
physical viscosity is larger than the unavoidable numerical viscosity. 



3.4 Radiative Transfer 

In many astrophysical problems, like the reionization of the Universe or in first star for- 
mation, one would like to self-consistently follow the coupled system of hydrodynamical 
and radiative transfer equations, allowing for a proper treatment of the backreaction of 
the radiation on the fluid, and vice versa. Unfortunately, the radiative transfer equation 
is in itself a high-dimensional partial differential equation that is extremely challenging 
to solve in its full generality. The task is to follow the evolution of the specific intensity 
field, 

- 3^ + n • V/ v = -K v / V + Jv , (56) 
c at 

which is not only a function of spatial coordinates, but also of direction n and frequency 
of v. Here j v is a source function, and K v an absorption coefficient, which provide an 
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implicit coupling to the hydrodynamics. Despite the formidable challenge to develop ef- 
ficient numerical schemes for the approximate solution of the radiation-hydrodynamical 
set of equations, the development of such codes based on SPH-based techniques has 
really flourished in recent years, where several new schemes have been proposed. 

Perhaps the simplest and also oldest approa ches are based on flux-lim ited diffusion 
approximations to radiative transfer in SPH dWhitehouse & Batel 120041: Whitehouse, 



tl [-/[-/i. W/Y All AH. LAW IAD IV/ A LlVl ILll 1 V \^ 11 UlllllVl ^^A^^^^^^^LJ^^^^ULA LUllUtMU V-*-- t-J LI l , ^\J\J^ * IT ill 

Bate & Monaghan.l2005l:IWhitehouse & BateLl2006l :lFrver. Rockefeller & WarrenL 120061: 



Viau. Basti en & Cha, 2006; Forsan et al.,|2009|). A refinement of this strategy has been 
proposed by IPetkova & Springel (2009), based on the optically thin variable tensor ap- 
proximations of iGnedin & Abell (1200 lb . In this approach, the radiative transfer equation 
is simplified in moment-based form, allowing the radiative transfer to be described in 
terms of the local energy density J v of the radiation, 

J v = ^-[l v da. (57) 



471 

This radiation energy density is then transported with an anisotropic diffusion equation 

1 dJ v d ( 1 dJ v h ij 



c 



dt drj \K V d 



Ky/v+jv (58) 



where the matrix h'i is the so-called Eddington tensor, which is symmetric and normal- 
ized to unit trace. The Eddington tensor encodes information about the angle depen- 
dence of the local radiation field, and in which direction, if any, it prefers to diffuse 
in cas e the local radiation field is highly anisotropic. In the ansatz of IGnedin & Abell 
d200lh . the Eddington tensors are simply estimated in an optically thin approximation as 
a 1/r 2 -weighted sum over all sources, which can be calculated efficiently for an arbitrary 
numb er of sources with tech niques familiar from the calculation of gravitational fields. 
Petkova & Springel (2009) derive an SPH discretization of the transport part of this 



radiative transfer approximation (described by the first term on the right hand side) as 

^=£*y(A>-Ai), (59) 
j 



where the factors wn are given by 



2cm ij rj j h i jV i W i j 

— 5 ' ( 6 °) 



and TV,- = c J v / (hp lanck v ) is the photon number associated with particle i. The tensor 

& = ^h-iTr(h) (61) 

is a modified Eddington tensor such that the SPH discretization in Eqn. (l60l corresponds 
to the correct anisotropic diffusion operator. The exchange described by Eqn. d59l is 
photon conserving, and since the matrix Wjj is symmetric and positive definite, rea- 
sonably fast iterative conjugate gradient solvers can be used to integrate the diffusion 
problem implicitly in time in a stable fashion. 

A compl etely different app r oach t o combine radiative transfer with SPH has been de- 
scribed by IPawlik & Schavel d2008h . who transport the radiation in terms of emission 
and reception cones, from particle to particle. This effectively corresponds to a coarse 
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discretization of the solid angle around each particle. Yet another scheme is given by 



Nayakshin. Cha & Hobbsl ((2009), who implemented a Monte-Carlo implementation of 



photon transport in SPH. Here the photons are directly implemented as virtual parti- 
cles, thereby providing a Monte-Carlo sampling of the transport equation. Although this 
method suffers from the usual Monte-Carlo noise limitations, only few simplifying as- 
sumptions in the radiative transfer equation itself have to be made, and it can therefore 
be made i ncreasingly more accurate simply by using more Monte-Carlo photons. 

Finally. lAltay. Croft & Pelupessyl ( 20081 suggested a ray-tracing scheme within SPH, 
which essentially implements ideas that are known from so-called long- and short char- 
acteristics methods for radiative transfer around point sources in mesh codes. These 
approaches have the advantage of being highly accurate, but their efficiency rapidly de- 
clines with an increasing number of sources. 



3.5 Relativistic Dynamics 



It is also possible to derive SPH eq uations for relativistic dynamics from a variational 
principle (IMonaghan & Price! 1200 lb . both for special and general relativistic dynamics. 



This is more elegant than alternative derivations a nd avoids some pr oblems inherent in 
other approaches to relativistic dynamics in SPH (Monaghan. 1992; Laguna, Miller & 
Zurek. ll993l:ISiegler & RiffertL l200d: iFaber & Rasiol l2000l : lAval etail l200lh . 
The variational method requires discretizing the Lagrangian 



T^U v U v dV, 



(62) 



where U^ 1 is the four- velocity, and 7^ v = (P + e)U v U v +-Pri jUV is the energy momentum 
tensor of a perfect fluid with pressure P and rest-frame energy density 

e = n moc 2 ( 1 + - 



(63) 

Here a (—1,1,1,1) signature for the flat metric tensor Tj^ is used, n is the rest frame 
number density of particles of mean molecular weight mo, and u is the rest-frame thermal 
energy per unit ma ss. In the following, we set c = jtiq = 1. 

Rosswog (2009) gives a detailed derivation of the resulting equations of motion when 
the discretization is done in terms of fluid parcels with a constant number v, of baryons 
in SPH-particle i and when variable smoothing lengths are consistently included. In the 
special-relativistic case, he derives the equations of motion as 



ds,- 
~d7 



fi^ViWM+fj^vwjihj) 



where the generalized momentum s,- of particle i is given by 

s; = ym (l + Ui + - 1 j , 



(64) 



(65) 



and y is the particle's Lorentz factor. The baryon number densities Nj in the comput- 
ing frame are estimated just as in Eqn. @, except for the replacements p, — > Ni and 



nij — > V;. Similarly, the correction factors /; are defined as in Eqn. (1231 ) with the same 
replacements. 

One also needs to complement the equation of motion with an energy equation of the 
form 



de,- x 

dT = -^ 



i j 



(66) 
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where 

1 + Uj 

e,- = v r si+ (67) 

Yi 

is a suitably defined relativistic energy variable. These equations conserve the total 
canonical momentum as well as the angular momentum. A slight technical complica- 
tion arises for this choice of variables due to the need to recover the primitive fluid vari- 
ables after each timestep from the updated integration variables s, and e,, which requires 
finding the root of a non-linear equation. 



4 APPLICATIONS OF SPH 



The versatility and simplicity of SPH has led to a wide range of applications in astron- 
omy, essentially in every field where theoretical research with hydrodynamical simula- 
tions is carried out. We here provide a brief, necessarily incomplete overview of some 
of the most prominent topics that have been studied with SPH. 



4.1 Cosmological Structure Formation 



Among the most important successes of SPH in cosmology are simulations that have 
clarifi ed the origin of the Lyman-q forest in the absorption spectra to distant quasars 
(e.g. iHernquist et all 119961 : iDave et all Il999h . which have been instrumental for test- 
ing and interpreting the cold dark matter cosmology. Modern versions of such calcula- 
tions use detailed chem i cal mo dels to study the enrichment of the intragalactic medium 
(lOppenheimer & Pavel I2006T) . and the nature of the so-c alled warm-hot IG M, which 
presumably contains a large fraction of all cosmic baryons (IDave et alll200ll) . 

Cosmological simulations with hydrodynamics are also used to study galaxy forma- 
tion in its full cosmological context, directly from the primordial fluctuation spectrum 
predicted by inflationary theory. This requires the inclusion of radiative cooling, and 
sub-resolution models for the treatment of star formation and associated energy feed- 
back processes from supernovae explosions or galactic wind s. A large variety of such 
models have been proposed (e.g. ISpringel & HernquistL l2003l) . some also involve rather 
substantial changes of SPH, for example the 'decoupled' version of SPH bv Marri & 
White (|2003h for the treatment of multiphase structure in the ISM. Despite the uncer- 
tainties such modelling involves and the huge challenge to numerical resolution this 
problem entails, important theoretical results on the clustering of galaxies, the evolution 
of the cosmic star formation rate, or the efficiency of galaxy formation as a function of 
halo mass have been reached. 

One important goal of such simulations has been the formation of disk galaxies in a 
cosmological context, which has however proven to represent a significant challenge. In 
recent years, substantial progress in this area has been achieved, however, where some 
simulated galaxie s have become quite close now to the morphology of real spiral galaxies 
(IGovernato et all 120071 : IScannapieco et all l2008h . Still, the ratio of the stellar bulge to 
disk components found in simulations is typically much higher than in observations, and 
the physical or numerical origin for this discrepancy is debated. 

The situation is better for clusters of galaxies, where the most recent SPH simulations 
have been quite successful in reproducing the primary scal i ng relations that are ob served 
dBorgani et all 120041 : Ipuchwein. Siiacki & Springe! 120081: iMcCarthv et alll2O10h . This 
comes after a long struggle in the literature with the cooling-flow problem, namely the 
fact that simulated clusters tend to efficiently cool out much of their gas at the centre, an 
effect that is not present in observations in the predicted strength. The modern solution 
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is to attribute this to a non-gravitational heat source in the centre of galaxies, which 
is thought to arise from active galactic nuclei. This important feedback channel h as 



recently been incorporated in SPH simulations of galaxy clusters (ISijacki et al. , 2008). 



SPH simulations of cosmic large-scale structu re and galaxy cluste rs are also regu 



larly used to stu dy the Sunyaev-Zeldovich effect ( da Silva et all l2000h, the bu ild-up of 
magnetic fields (|Dolag. Bartelmann & Leschlll999ll2002!:lDolag et all l2005a|) . or even 
the production of non-thermal particle populations in the form of cosmic rays (Pfrom- 
mer et al., 120071 : lJubelgas et all l2008|) . The latter also required the development of a 
technique to detect sho ck waves on the fly in SPH, and to measure their Mach number 
dPfrommer et alll2006h . 



4.2 Galaxy Mergers 

The hierarchical theory of galaxy formation predicts frequent mergers of galaxies, lead- 
ing to the build-up of ever larger systems. Such galaxy interactions are also promi- 
nently observed in man y systems, such as the 'A ntennae Galaxies' NGC 4038/4039. 
The merger hypothesis ( Toomre & Toomrel 1972|) suggests that the coalescence of two 
spiral galaxies leads to an elliptical remnant galaxy, thereby playing a central role in ex- 
plaining Hubble's tuning fork diagram for the morphology of galaxies. SPH simulations 
of merging galaxies have been instrumental i n understanding this process . 

In pioneering work bv Hernquist (1989), Barnes & Hernquist (1991) and Mihos & 
Hernquist (119941 1 19961) . the occurrence of central nuclear starbursts during mergers was 
studied in detail. Recently, the growth of supermassive b lack holes has been added to 
the simulations (ISpringel. Di Matteo & Hernquist , 120051) . making it possible to study 
the co-evolution of black holes and stellar bulges in galaxies. Di Matteo, Springel & 
Hernquist (2005) demonstrated that the energy output associated with accretion regu- 
lates the black hole growth, and establishes the tight observed relation between black 
hole masses and bulge velocity dispersions. Subsequently, these simulations have been 
used to develop comprehensive models of spheroid for mation and for interpreting t he 
evolution and properties of the cosmic quasar population (IHopkins et all |2005[ 12006). 



4.3 Star Formation and Stellar Encounters 

On smaller scales, many SPH-based simulations have studi ed the fragm entation of molec- 
ular clouds, star formation, and the initial mass function (Bate, 1998; Klessen, Burkert 
& Bate. ll998l:lBate & BonnellLl2005l : ISmith. Clark & BonnellL bOOglTThis includes also 
simulations of the formatio n of the first stars in the Un i verse (IBromm. Coppi & Larson , 
20021: lYoshida et ail 120061: IClark. Glover & Klessenl 120081) . which are thought to be 
very special objects. 

In this area, importa nt numerical issues related to SPH have been examined as well. 
Bate & Burkertl (119971) show that the Jeans mass needs to be resolved to gua rantee reli- 
able re sults for gravitational fragmentation, a result that was strengthened by Iwhitworthl 
(119981) . who showed that provided h ~ £ and the Jeans mass is resolved, only physical 
fragmentation should occur. Special on-the-fly particle splitting methods have been de- 
veloped in order to guarantee that t he Jeans ma ss stays resolved over the coarse of a 
simulation (|Kitsionas & Whitworthl l200l 120071). Also, the ability of SPH to repre sent 
driven isothermal turbulence has been studied (IKlessen. Heitsch & Mac Lowl 1200 0). 

Pion eeri ng work on stella r collisions with SPH has been performed by iBenz & Hills 
dl987h and Benz et al. (1990). Some of the most recent work in this area studied quite so- 
phisticated problems, such as relativistic neutron star mergers with an approximate treat- 
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ment of general relativity and sophisticated nuclear matter equations of state (Oechslin, 
Janka & Marek, 120071) . or the triggering of su b-luminous supernovae type la through 
collisions of white dwarfs ( Pakmor et allbOlOT) . 



4.4 Planet Formation and Accretion Disks 

SPH has been successfully emplo yed to study giant planet formation through the frag 
mentation of protoplanetary disks (|Mayer et all I2002L I2004Q . and to ex plore the interac 



tion of high- and lo w-mass planets embedded in protoplanetary discs (IB ate et all 12003 
Lufkinetal!l2004h . 



1995), 



1997), 



A number of works have also used SPH to m odel accretion disks (Simpsonl 
including studies of spiral shocks in 3D disks (I Yukawa. Boffin & MatsudaL 
accretion o f gas onto b lack hol es in viscous accretion disks (Lanzaf ame, Molteni & 
Chakrabarti, 1998), or cataclysmic variable systems dWoodetalll2005h . 

Another interesting application of SPH lies in simulations of materials with differen t 
equations of state, corresponding to rocky and icy materials ( Benz & Asphaug . 1999|) . 
This allows the prediction of the outcome of collisions of bodies with sizes from the 
scale of centimeters to hundreds of kilometers, which has important implications for the 
fate of small fragments in the Solar system or in protoplanetary disks. Finally, such 
simulation techniques allow calculations of the collision of protoplanets, culminating 
in numerical simulations that showed how the Moo n might have formed by an impac t 
between protoearth and an object a tenth of its mass (IBenz. Slattery & Cameron. Il986h . 



5 CONVERGENCE, CONSISTENCY AND STABILITY OF SPH 

There have bee n a few code comparisons in the literature between SPH and Eulerian 
hydrodynamics ( Frenk et all 1 1 9991 : lo ' Shea et al. . 2005 : Tasker et al. . 2008 ). but very few 
formal studies of the accuracy of SPH have been carried out. Even most code papers on 
SPH report only circumstantial evidence for SPH's accuracy. What is especially missing 
are rigorous studies of the convergence rate of SPH towards known analytic solutions, 
which is ultimately one of the most sensitive tests of the accuracy of a numerical method. 
For example, this may involve measuring the error in the result of an SPH calculation 
in terms of an LI error norm relative to a known solution, as it is often done in tests 
of Eulerian hydrodynamics codes. There is no a priori reason why SPH should not 
be subjected to equally sensitive tests to establish whether the error becomes smaller 
with increasing resolution (convergence), and whether the convergence occurs towards 
the correct physical solution (consistency). Such tests can also provide a good basis to 
compare the efficiency of different numerical approaches for a particular type of problem 
with each other. 

In this section, we discuss a number of tests of SPH in one and two dimensions in 
order to provide a basic characterization of the accuracy of 'standard SPH' as outlined 
in Section [2] Our discussion is in particular meant to critically address some of the 
weaknesses of standard SPH, both to clarify the origin of the inaccuracies and to guide 
the ongoing search for improvements in the approach. We note that there is already a 
large body of literature with suggestions for improvements of 'standard SPH', ranging 
from minor modifications, say in the parametrization of the artificial viscosity, to more 
radical changes, such as outfitting SPH with a Riemann solver or replacing the kernel- 
interpolation technique with a density estimate based on a Voronoi tessellation. We shall 
summarize some of these attempts in Section [6l but refrain from testing them in detail 
here. 
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Figure 1: Convergence rate for a traveling sound wave in ID, calculated in SPH without artificial 
viscosity. The solution (symbols) converges with second order accuracy as LI °c N~ 2 (this power 
law is shown with dashed lines), for different numbers of neighbours. 



We note that for all test results reported below, it is well possible that s mall modi- 
ficatio ns in the numerical parameters of the code that was used (GADGET2, ISpringelL 
20051) may lead to slightly improved results. However, we expect that this may only 
reduce the error by a constant factor as a function of resolution, but is unlikely to sig- 
nificantly improve the order or convergence. The former is very helpful of course, but 
ultimately represents only cosmetic improvements of the results. What is fundamentally 
much more important is the latter, the order of convergence of a numerical scheme. 



5.1 One-dimensional Sound Waves 



Based on analytic reasoning, Rasiol ( 2000 ) argued that convergence of SPH for one- 
dimensional sound waves requires increasing both the number of smoothing neighbours 
and the number of particles, with the latter increasing faster such that the smoothing 
length and hence the smallest resolved scale decreases. In fact, he showed that in this 
limit the dispersion relation of sound waves in ID is correctly reproduced by SPH for all 
wavelengths. He also pointed out that if the number of neighbours is kept constant and 
only the number of particles is increased, then one may converge to an incorrect physical 
limit, implying that SPH is inconsistent in this case. 

We note that in practice the requirement of consistency may not be crucially important, 
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Pi 


Pi 


Vl 


P2 


Pi 


V2 


Problem 1 


1.0 


1.0 


0.0 


0.125 


0.1 


0.0 


Problem 2 


1.0 


0.4 


-2.0 


1.0 


0.4 


2.0 


Problem 3 


1.0 


1000.0 


0.0 


1.0 


0.01 


0.0 



Table 1: Parameters of the one-dimensional Riemann problems (for y= 1 .4) examined here. 



provided the numerical result converges to a solution that is 'close enough' to the correct 
physical solution. For example, a common error in SPH for a low number of neighbours 
is that the density estimate carries a small bias of up to a percent or so. Although this 
automatically means that the exact physical value of the density is systematically missed 
by a small amount, in the majority of astrophysical applications of SPH the resulting 
error will be subdominant compared to other errors or approximations in the physical 
modelling, and is therefore not really of concern. In our convergence tests below we 
therefore stick with the common practice of keeping the number of smoothing neigh- 
bours constant. Note that because we compare to known analytic solutions, possible 
errors in the dispersion relation will anyway be picked up by the error measures. 

We begin with the elementary test of a simple acoustic wave that travels through a 
periodic box. To avoid any wave steepening, we consider a very small wave amplitude 
of Ap/p = 10~ 6 . The pressure of the gas is set to P = 3/5 at unit density, such that 
the adiabatic sound speed is c s = 1 for a gas with y= 5/3. We let the wave travel once 
through the box of unit length, and compare the velocity profile of the final result at time 
t = 1.0 with the initial conditions in terms of an LI error norm. We define the latter as 

L1 = ^Z>«— v (*«-)l, ( 6§ ) 

where N is the number of SPH particles, v, is the numerical solution for the velocity 
of particle i, and v(jc,-) is the expected analytic solution for the problem, which is here 
identical to the initial conditions. Note that we deliberately pick the velocity as error 
measure, because the density estimate can be subject to a bias that would then completely 
dominate the error norm. The velocity information on the other hand will tell us more 
faithfully how the wave has moved, and whether it properly returned to its original state. 

In Figure [Q we show results for the LI error norm as a function of the number N 
of equal-mass particles used to sample the domain, for different numbers of smoothing 
neighbours. Interestingly, the results show second-order convergence of the code, with 
LI oc N~ 2 , as expected in a second-order accurate scheme for smooth solutions without 
discontinuities. This convergence rate is the same as the one obtai ned for this pr oblem 
with standard state-of-the-art second-order Eulerian methods (e.e. ISr)ringelL[2 009: Stone 



et al.. 120081) . which is reassuring. Note that in this test a larger number of neighbours 
simply reduces the effective spatial resolution but does not lead to an advantage in the 
LI velocity norm. Nevertheless, the density is more accurately reproduced for a larger 
number of neighbors. Because the sound speed depends only on temperature (which is 
set as part of the initial conditions), the density error apparently does not appreciably 
affect the travel speed of acoustic waves in this test. 

5.2 One-dimensional Riemann Problems 



Let us now consider a few Riemann problems in one dimension. Their initial conditions 
are characterized by two piece-wise constant states that meet discontinuously at x = 0.5 
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Figure 2: Different one-dimensional Riemann problems, calculated with a resolution of 100 
points in the unit domain, and 7 smoothing neighbours. The three columns show results for the 
initial conditions of the problems 1, 2 and 3 as specified in Table Q] Symbols represent the SPH 
particles, solid lines the analytic solutions for density, velocity and pressure, from top to bottom. 



at time t = 0. The subsequent evolution then gives rise to a set of self-similar waves, 
containing always one contact wave, which is sandwiched on the left and the right by 
either a shock wave or a rarefaction wave. The ability to correctly reproduce the non- 
linear outcome of arbitrary Riemann problems is the backbone of any hydrodynamical 
method. 

In Figure [2] we show SPH results for three different Riemann-type problems, with 
initial conditions characterized by triples of density, pressure and v elocity for e ither side, 
as listed in Table Q] For all these problems (which are taken from lToroUl997h . y= 1.4 
has been assumed. Problem 1 gives rise to a comparatively weak shock, which has been 
studied in very similar fo rm (but with a different value of P? = 0.1795) in a number 
of previous test s of SPH dHernauist & KatzL Il989l: iRasio & ShapiroL Il99ll: Wadslev. 
Stadel & Ouinn, l2004l : ISpringelLl2005[ among others). In the second problem, the gas 
is suddenly ripped apart with large supersonic velocity, giving rise to a pair of strong 
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Figure 3: Convergence rate of one-dimensional Riemann problems (as labeled), calculated with 
a resolution of 100 points in the unit domain, using 7 smoothing neighbours. The dashed lines 
indicate aLl« A^ 1 scaling of the error. 



rarefaction waves. Finally, problem 3 is a again a Sod shock type problem but involves 
an extremely strong shock. 

The results shown in Figure [2] are based on 100 points with initially equal spacing in 
the unit domain, 7 smoothing neighbours, and a standard artificial viscosity setting. All 
problems are treated qualitatively correctly by SPH, with some inaccuracies at the con- 
tact discontinuities. Characteristically, shocks and contact discontinuities are broadened 
over 2-3 smoothing lengths, and rarefaction waves show a small over- and underestimate 
at their high- and low-density sides, respectively. Also, there is a pressure "blip" seen at 
the contact discontinuity. However, the properties of the postshock flow are correct, and 
the artificial viscosity has successfully suppressed all postshock oscillations. Also, the 
errors become progressively smaller as the resolution is increased. This is seen explicitly 
in Figure [3l where we show results for the LI error norm for the velocity as a function of 
resolution. The error declines as LI <=c TV -1 , which is expected due to the reduced order 
of the scheme at the discontinuities. The same convergence rate for the shock-problem is 
obtained with Eulerian approaches, as they too exhibit only first order accuracy around 
discontinuities in the solutions. 

Often the discussion of the numerical accuracy of shocks focuses on the sharpness 
with which they are represented, and on that basis SPH has frequently been portrayed 
as being inferior in comparison with Eulerian methods. However, it should be noted 
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Figure 4: Velocity profile and convergence rate of a two-dimensional strong shock problem. In 
the left panel, symbols show the velocities of all SPH particles, compared to the analytic result 
(solid line). In the right panel, the open circles show the LI error for a binned result (using N 
bins in the ^-direction in each case), the filled circles give the error when individual particles are 
considered. In this 2D problem, the ideal A^ 1 convergence rate (dotted lines) of ID simulations 
is reduced to approximately LI °< N~ QJ (dashed lines). 

that in both approaches the numerical shock width is always many orders of magnitude 
larger than the true width of the physical shock layer. What is much more relevant than 
the width is therefore that the properties of the post-shock flow are correct, which is the 
case in SPH. Also, the numerical thickness of the shocks may be reduced arbitrarily by 
using more particles, and is hence entirely a matter of resolution. We also point out that 
contact discontinuities can often become quite broad in Eulerian methods; in fact, their 
width increases with growing advection speed, unlike in the Galilean-invariant SPH. As 
far as one-dimensional hydrodynamics is concerned, SPH can hence be characterized 
as being quite accurate, and in particular, its convergence rate appears competitive with 
second-order Eulerian schemes. 

5.3 Two-dimensional Shock Waves 

Multi-dimensional hydrodynamics adds much additional complexity, such as shear flows, 
fluid instabilities, and turbulence. We here first briefly examine whether a problem with 
one dimensional symmetry, again a Riemann problem, can be equally well represented 
with SPH in multiple dimensions. For definiteness, we study again the strong shock of 
problem 3 from the previous section, but this time with the j-dimension added to the 
initial set-up. This is taken as a regular Cartesian N xN particle grid in the unit square, 
with periodic boundaries in the ^-direction, and reflective boundaries in the x-direction. 

In the left panel panel of Figure 01 we show the x velocities of all particles of a run at 
100 x 100 resolution and compare them to the analytic solution. Relative to the corre- 
sponding ID-result shown in the middle-right panel of Figure [2l the primary difference 
is a considerably increased 'noise' in the particle velocities. This noise is a generic fea- 
ture of multi-dimensional flows simulated with SPH. It can be reduced by an enlarged 
artificial viscosity or a larger number of neighbours, but typically tends to be much larger 
than in ID calculations. Presumably, this is just a consequence of the higher degree of 
freedom in the particle motion. 
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Unsurprisingly, the noise has a negative impact on the convergence rate of planar 
shocks in multi-dimensional SPH. This is shown for the same problem in FigureHl where 
we measure LI oc N~ QJ instead of the ideal LI oc N~ L0 . It also does not appear to help 
to first bin the particles and then to compute the LI error as a difference between the 
averaged and the analytical result. Although this reduces the absolute size of the error 
substantially and is certainly warranted to eliminate the intrinsic noise, it does not affect 
the convergence rate itself (see Fig. HI). Nevertheless, the noise does not destroy the 
principal correctness of the solution obtained with SPH for multi-dimensional planar 
Riemann problems. 



5.4 Shear Flows 

We now a discuss a yet more demanding two-dimensiona l problem, the vortex test of 



Gresho (Gres ho & Chan, Il990l : iLiska & W endroffl l2003h . It consists of a triangular 



azimuthal velocity profile 

0< r <0.2 

v^{r) = { 2-5r for 0.2 < r < 0.4 (69) 
r>0A 

in a gas of constant density equal to p = 1 and adiabatic index of y = 5/3. By adopting 
the pressure profile 




P{r) 



' 5 + 25/2r 2 for < r < 0.2 

(70) 



9 + 25/2r 2 



20r + 41n(r/0.2) for 0.2 < r < 0.4 
3 + 41n2 for r > 0.4 



the centrifugal force is balanced by the pressure gradient, and the vortex becomes inde- 
pendent of time. 

In Figure [51 we compare the results for the azimuthal velocity profiles at time t = 1 .0 
for three different runs, carried out with 80 x 80 particles in the unit domain for different 
settings of the artificial viscosity. The left-most panel shows the outcome for a standard 
viscosity of a = 1 .0, the middle panel is for a = 0.05 and the right panel did not use any 
artificial viscosity at all. We see that in all three cases substantial noise in the velocity 
profile develops, but it is clearly largest in the simulation without viscosity. On the other 
hand, one can see that the average velocity profile of the run without viscosity is actually 
closest to the expected stationary solution (blue lines), while the standard viscosity run 
already shows a reduced angular frequency in the solid-body part of the rotation in the 
inner part of the vortex. Despite the use of the Balsara switch in this problem, the 
velocity noise, and as a consequence the noisy estimates of divergence and curl, have 
produced enough residual viscosity to lead to appreciable angular momentum transport. 

This is corroborated by the results of convergence tests for this problem. In Fig- 
ure [6j we consider the LI error norm of the binned azimuthal velocity profile (based on 
N/2 bins, linearly placed in the radial range to 0.5). Interestingly, the two runs with 
non-vanishing viscosity do not converge to the analytic solution with higher resolution. 
Instead, they converge to a different solution. If we approximately identify this solution 
as the highest resolution result in each case, then we see that the lower resolution calcu- 
lations converge to it with LI oc /y~° 7 . What is important to note is that the effects of 
the artificial viscosity on the shear flow do not diminish with higher resolution. Instead, 
the solutions behave as if one had simulated a fluid with some residual shear viscosity 
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Figure 5: Velocity profile in the Gresho vortex test at time t = 1.0, at a base resolution of 
80 x 80, for different viscosity settings. Here 'low viscosity' corresponds to a = 0.05, and 
'normal viscosity' to a = 1.0, with the Balsara switch enabled in both cases. The small dots are 
the azimuthal velocities of individual SPH particles, while the red circles show binned results. 
The blue triangular shaped profile is the analytic (stationary) solution. 



instead. Only the run without viscosity appears to converge to the correct physical so- 
lution, albeit at a very low rate for high N. In any case, it appears that the convergence 
rate of SPH for this problem is considera bly worse than t he LI oc N~ lA measured for a 
moving-mesh code and Eulerian codes by lSpringel d2009h . 



5.5 Contact Discontinuities and Fluid Instabilities 

An issue that has created considerable attention recently is the question of whether SPH 
can properly resolve fluid instabilities, such as the Kelvin-Helmholtz (KH) instability in 
shear flows. For a set-up with equal particle masse s and a sharp initial density contrast of 
1 : 2 at a contact discontinuity, lAgertz et all (|2007h did not observe any gro wth of the in 



stabili ty, whereas for a vanishing density jump the fluids would start to mix. lAgertz et al. 



(2007) attributed these problems to substantial errors in the pressure gradient estimates 
at the contact discontinuity. A number of recent studies have addressed this problem as 
well, prop o sing m odifications of SPH designed to improve the KH results. In partic- 
ular, Price ( 2008|) has suggested adding an artificial heat conduction to smooth out the 
ph ase boundary, which indeed impr oved the results. A different approach was followed 
by iRead. Hayfield & Agertz d2O10h . who employed a different kernel, a much enlarged 
number of neighbors, and a modified density estimation formula to obtain a better rep- 
resentation of mixing of different phases in shear layers. We here do not examine these 
modifications but want to shine more light onto the nature of the inaccuracies of standard 
SPH for the KH problem. 

We adopt the KH test parameters of Springe 1 (2009), i.e. a central phase in the region 
\y — 0.5 1 < 0.25 is given density P2 = 2 and velocity v x = 0.5, while the rest of the gas has 
density pi = 1 and velocity v x = —0.5, all at the same pressure of P = 2.5 withY= 5/3 . 
In order to soften the transition at the two interfaces, we follow Robertson et al. (2009) 
and impose a smooth transition in the vertical profiles of density p(y), shear velocity 
v x (y), and specific entropy A(y) = P/p(j) Y . For example, the initial density structure is 
adopted as 



p(y) = pi + 



P2-P1 



[1 + exp(-2(j - 0.25)/o)] [1 + exp(2(j - 0.75)/o)] ' 



(71) 



with a = 0.025, and similar adoptions are made for the other profiles. The transition 
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Figure 6: Convergence rate of the Gresho problem, in terms of the binned azimuthal velocity, 
for different viscosity settings. Here the open symbols measure the LI error towards the highest 
resolution result, while the filled symbols measure the difference with respect to the analytic 

convergence. 



solution. The dotted line indicates a LI °< N 7 



removes the sharp discontinuities from the initial conditions, such that they can in prin- 
ciple be fully resolved by the numerical scheme. This regularization makes the problem 
well-posed for convergence studies, which is particularly important for the KH instabil- 
ity, where smallest scales grow fastest. 

For simplicity, we consider a Cartesian particle grid in the unit domain with periodic 
boundaries, and trigger the instability by applying the velocity perturbation throughout 
the y-domain, as 

v y (x) = vq sin(fcx), (72) 

with an initial amplitude of vo = 0.01 and a wave number of k = 2 x (2n) /L. We can 
easily study the growth of this velocity mode by measuring its amplitude through a 
Fourier transform of the two-dimensional v y velocity field. 

In Figure |7J we show density maps of the evolved fields at time t = 2.0 for different 
SPH simulations, using 27 neighbours and no artificial viscosity. As can be seen, the 
KH instability clearly develops and produces the characteristic wave-like features. The 
density field is visibly noisy though, as a result of the absence of artificial viscosity. If 
the latter is added, the noise is substantially reduced, but the KH billows look much more 



30 



V. Springel 




anemic and have a reduced amplitude. 

An interesting question is whether the KH instability actually grows with the right 
rate in these SPH calculations. This is examined in Figure [H where the amplitude of 
the excited v y mode of the velocity field is shown as a function of time, for the case 
without viscosity. The growth is compared to the expected exponential growth rate v y <x 




Figure 8: Growth rate of KH instabilities in SPH test simulations, for different resolutions as 
labeled. The left panel gives the amplitude of the initially excited velocity mode, compared to 
the expected linear growth rate (dashed line). The right panel shows the velocity dispersion of 
the SPH particles in the x-direction as a function of time, for a thin layer close to the midplane 
of the box. This dispersion is an indication of the level of velocity noise in the calculation. 
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Figure 9: The same as in Figure [8] but for simulations that include an artificial viscosity of 
moderate strength with a = 0.25. 



exp(?/iKH) shown as a solid line, where 

t kh = 1 1 : . (73) 

|v 2 — vi| A:- v /pip 2 

is the KH growth timescale for an inviscid gas. After an initial transient phase (which is 
expected because the initial perturbation was not set-up self-consistently), the two low- 
resolution calculations with 32 2 and 64 2 particles do actually reproduce this growth rate 
quite accurately for some time, but then the growth rate suddenly slows down consider- 
ably. Curiously, the higher resolution simulations never quite reach the expected growth 
rate. So what is going on here? 

A hint is obtained by the results shown in the right panel of Fig.[8j where the v x veloc- 
ity dispersion of the particles in a narrow strip around y ~ 0.5 is shown, as a function of 
time. We in principle expect this to remain close to zero until very late in the evolution. 
However, what we actually observe is a sudden increase of this dispersion, when the reg- 
ularity and coherence of the initial Cartesian particle grid is finally lost. At this moment, 
SPH develops its velocity noise that is characteristic even of smooth flows. Interestingly, 
we observe that the development of this velocity noise closely coincides in time with the 
termination of the correct linear growth rate. The natural interpretation is that as long as 
the particles are still quite regularly ordered, the analytic linear KH growth rate is calcu- 
lated accurately by SPH, but at late times the developing noise introduces a substantial 
degradation of the accuracy that reduces the proper growth rate. 

This suggests that perhaps the use of the right amount of artificial viscosity may yield 
an accurate growth rate even at late times if it prevents the detrimental noise effects. 
In Figure |9l the equivalent results are shown with artificial viscosity included. Here 
the behaviour is rather similar initially. As long as the regular order of the particles 
is still approximately maintained, the artificial viscosity is not changing anything as it 
is reduced by the Balsara switch to a very low level. Once the initial order decays, 
the velocity dispersion suddenly increases, first to the same level that one would have 
without viscosity. Only after some time, the viscosity damps out this noise, but this 
comes at the price of also damping the growth rate even further, since the fluid behaves 
now in a slightly viscous fashion. 

We conclude that whereas KH fluid instabilities do occur in SPH, calculating them 
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with high accuracy is a challenge due to the significant noise present the multi-dimensional 
velocity field. Getting rid off this noise with artificial viscosity tends to reduce the KH 
growth rate below the analytic expectation, making it hard to achieve truly inviscid be- 
havior. More work on this problem is therefore clearly warranted in the future. 



5.6 Surface Tension 

In the standard energy and entropy conserving formulation of SPH, two phases of parti- 
cles with different specific entropies tend to avoid mixing, simply because this is energet- 
ically disfavored. To understand the origin of this effect, let us consider a simple virtual 
experiment where a box of volume V is filled in one half with gas of density p i , and in 
the other half with gas of density P2. We assume an equal pressure P, such that the SPH 
particles in the different phases will have specific entropies Ay = P /p\ and A2 = P/ pi- 
Then the total thermal energy in the box will be .Etherm = u\M\-\- U2M2 = PV /(y — 1). 
The numerical estimate for a particular SPH realization may be slightly different from 
this due to intermediate density values in the transition layer, but this effect can be made 
arbitrarily small if a large number of particles is used. Now imagine that we rearrange 
the particles by homogeneously spreading them throughout the volume, in some regu- 
lar fashion (say as two interleaved grids), but keeping their initial entropies. Provided 
the SPH smoothing lengths are large enough, each particle will then estimate the mean 
density p = (pi +p2)/2 as its new density, throughout the volume. The new ther- 
mal energy per unit mass for particles of species 1 will then be u\ = Aip y ~ 1 /(y— 1) 
and similarly for particles of species 2. As a result, the new thermal energy becomes 

E Lrm = Wfr- 1)] [(pi + P2)/2] Y_1 [p^ + P^] A which is larger than the orig- 
inal energy if pi 7^ P2- Clearly then, the mixing of the particle set in this fashion is 
energetically forbidden, and an energy conserving code will resist it. This resistance 
appears as an artificial surface tension term in SPH. 

We note that the mixing can be accomplished at constant thermal energy, but this 
requires that the entropies of the particles and hence their temperature estimates are made 
equal. The final entropy is then A = P /p y , which corresponds to a total thermodynamic 
entropy that is larger than in the unmixed state. Mixing the phases in this irreversible 
fashion hence requires creating the relevant amount of mixing entropy, but for this no 
source is foreseen in the entropy- conserving formulation of standard SPH. Including 
artificial thermal conduction (as in Pried. 2008h is one interesting approach to address 



this problem, as this process equilibrates the temperatures while conserving energy and 
increasing the entropy. 

The presence of some level of s urface tension in SPH c an be demonstrated experimen- 



tally through simple settling tests (|HeB & SpringelL |201fJ). For example, one possibility 



is to set-up an overdense spherical region and let it relax to an equilibrium distribution 
(this can be done by keeping the particle entropies fixed, and by adding an artificial de- 
cay of the velocities). Then the pressure inside of the sphere in the final relaxed state is 
found to be slightly higher than outside, as it needs to be to balance the surface tension. 
In fact, according to the Young-Laplace equation, the expected pressure difference is 

AP =°(^ + ^)' <74) 

where a is the surface tension, and Ry and R2 are the two principal radii of curvature 
of the surface. In Figure [lOl we show the outcome of such an experiment in 2D, for 
densities P2 = 2 and pi = 1 realized with equal mass particles, a pressure of P = 2.5, 
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Figure 10: Equilibrium configurations of particles in 2D for two phases of different entropy, 
corresponding to a jump in density by a factor of 2. Both in the left and in the right case, the final 
pressure in the inside of the sphere is slightly higher than the one on the outside. This pressure 
difference offsets the spurious surface tension present at the interface. 



and 13 smoothing neighbors. Independent of whether the high or low density phase 
is arranged to be inside the sphere, the inner region shows a small pressure difference 
relative to the outer region, which is AP ~ 0.002 at an effective resolution of N 2 = 128 2 
for the lower density phase. For a radius of Ri = 0.4 of the sphere, this then implies a 
surface tension of a = 0.0025. We have confirmed that this surface tension varies with 
one-dimensional resolution as a°= 1 /N, i.e. it declines with higher spatial resolution. 
On the other hand, it increases for higher density ratios, and for a larger number of 
smoothing neighbours. 

One consequence of this tension is that SPH in principle may support capillary waves 
at interfaces, with dispersion relation go 2 = ok 3 / (pi + P2). In the context of the Kelvin- 
Helmholtz problem considered above, we note that surface tension may also modify the 
growth rate of small wavelength perturbations. In the presence of surface tension, the 
KH growth timescale becomes 



?kh — 

and waves with wavelength 



£ 2 pip2(v 2 -vi) 2 ok 3 



(pl+p2) 2 Pl+p2_ 



-1/2 

(75) 



k < A cr it = 271 -y (76) 

PlP2 (V2-Vl) Z 

will not grow at all. Using the value measured for a for the numerical set-up we consid- 
ered above, we obtain X C n t /d ~ 3.0/ (v2 — vi) 2 , where d is the mean particle separation. 
For a shear of V2 — vi = 1, we hence expect only waves with wavelengths up to a cou- 
ple of particle separations to be suppressed. In particular, this effect should not have 
influenced the KH tests carried out in the previous section. However, for small shear, 
the effect can be much more of a problem, and here may stabilize SPH against the KH 
instability. 
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Figure 11: Clumping instability in a 2D settling test, using 13 (left panel), 27 (middle panel), or 
55 (right panel) smoothing neighbours. In all three cases, the number of particles used (50 2 ) is 
the same. It is seen that for a large number of smoothing neighbours, particles clump into groups, 
thereby reducing the effective spatial resolution of the scheme. 



5.7 The Tensile Instability 

First described by ISwegld dl995h . the clumping or tensile instability is a well known 
nuisance in SPH that occurs if a large number of smoothing neighbours is used. In this 
case, it can happen that the net force between a close particle pair is not repulsive but 
attractive, simply because the kernel gradient tends to become shallow for close separa- 
tions such that the pair may be further compressed by other surrounding particles. As 
a result, particles may clump together, thereby reducing the effective spatial resolution 
available for the calculation. 

Settling tests in which an initially random distribution of particles with equal and 
fixed entropies is evolved towards an equilibrium state under the influence of a friction 
force are a quite sensitive tool to reveal the presence of this instability. We show an 
example of the outcome of such a test in Figure [TTJ where 50 2 particles were randomly 
placed in a periodic 2D box of unit length, and then evolved until a good equilibrium 
of the pressure forces was reached, corresponding to a low energy state of the system. 
The leftmost panel in Fig. [TT] shows the result for 13 smoothing neighbours, the middle 
panel for 27 neighbours, and the right panel for 55 neighbours. While for a low number 
of neighbours a nicely regular, hexagonal particle distribution is obtained, for a large 
number of smoothing neighbours many particle clumps are formed, strongly reducing 
the number of independent sampling points for the fluid. 

Taken at face value, this seems to make attempts to use a large number of neighbors 
to increase the accuracy of SPH calculations a futile exercise. However, we note that the 
clumping instability is essentially never observed in this extreme form in real applica- 
tions of SPH. This is because the dissipation added in the settling test is not present in 
real application, and there the dynamics is usually also not 'cold' enough to create many 
artificial SPH particle clumps or string-like configurations. Nevertheless, a number of 
suggestions hav e been made to fix the tensile instability through modifications of the 
SPH formalism. Steinmeta (199a) simply invoke an artificially steepened firs t derivative 
of the kernel to obtain a larger repulsive pressure force at small separations. iMonaghanl 
(2000) proposed instead an additional artificial pressure, which has the advantage of 
maintaining the conservativ e properties of SPH, but it introd uces some small error in the 
dispersion relation. Finally, IRead. Hayfield & Agertjd2010h suggest the use of a peaked 
kernel, thereby ensuring large repulsive pressure forces at small separations, at the price 
of a lower order of the density estimate. 
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6 ALTERNATIVE FORMULATIONS OF SPH AND FUTURE DIREC- 
TIONS 

There are many avenues for modifications of the standard SPH method outlined in Sec- 
tion |2j even though this usually implies giving up the elegance of the Lagrangian deriva- 
tion in favor of a more heuristic construction of the discretization scheme. However, 
what ultimately matters in the end is the achieved accuracy, and hence there is certainly 
a lot of merit in investigating such alternative schemes. We here briefly mention some 
of the proposed modifications. 



6.1 Modifications of the Kernel Estimation Scheme 



Changing the basic SPH kernel shape has been investigated in a number of studies (Mon- 



aghan, 19851: Fulk , 19961 : Cabezon. Garcia-Senz & Relafio . 2008), but generally no ker- 
nel has been found that performed significantly better than the cubic spline. In particular, 
higher order kernels may both be positive and negative, which invokes numerous stabil- 
ity problems. 

A few studies (jShapiro et al. , 19961 : lOwen et al. , 1998 ) proposed anisotropic kernels 
that adapt to local flow features. Together with the use of a tensorial artificial viscosity, 
this promises an optimum use of the available spatial resolution, and hence ultimately 
higher accuracy for a given number of particles. 

In the context of galaxy formation simulations and the so-called overcooling problem, 
it has been pointed out that SPH's density es timate in the vicinity of contact disconti- 
nuities may be problematic ( Pearce et al. . 1999}) . This led Ritchie & Thomas d200ll) to 
suggest a modified density estimate specifically designed to work better across contact 
discontinuities. To this extent, particles are weighted with their thermal energies in the 
density estimate, which works well across contacts. In particular, this removes the pres- 
sure blip that is often found in contact wa ves that develop in Riemann prob lems. This 
idea has also been recently reactivated by iRead. Hayfield & Agertz ( 201fj|) . combined 
with the suggestion of a much larger number of neighbors and a modified shape of the 
smoothing kernel. 

A more radical approach has been proposed by lBprve. Qmang & Trulser] ( 200 lb . who 
introduced the so-called regularized SPH formalism. In it the discretization errors of 
the ordinary kernel estimation are reduced by introducing a new concept of intermediate 
interpolation cells. They are also used for periodic redefinitions of the whole particle 
distribution to obtain a better numerical description of the dynamics, in particular of 
MHD related problems. 

Another attempt to reduce the errors associated with kernel s ums is given by Imaeda 
& Inutsuka (120021) in terms of their 'consistent velocity method'. The central idea of 
this scheme is to achieve consistency between the continuity equation and the density 
estimate. In some s ense this makes it reminiscent of the simpler idea of 'XSPH', orig- 
inally introd uced by Monaghan ( 1989|) and later derived in more accurate form from a 
Lagrangian (|Monag han. 2002). In this scheme the particles are moved with a SPH inter- 
polated version of the velocity, and not with the individual particle velocities. This has 
been used in some simulations to reduce problems of particle interpe netration. 



Finally, a higher order formulation of SPH has also been proposed (Oger et al., 2007) 



It has been shown to eliminate the tensile instability, but at the price of somewhat larger 
errors in the gradient estimates. 
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6.2 Improved Artificial Viscosity 

The idea to use a time- variable coefficient for the parametrization of the artificial viscos- 
ity is in general very interesting. It can help to restrict the viscosity to regions where it 
is really needed (i.e. in shocks), though it should be reduced to very small levels every- 
where else, ideally to the lowest level still consistent with keeping the velocity noise of 
the scheme under control. There is clearly still room for improvement of the switches 
used at present to increase or decrease the viscosity, as discussed in Section |23~ 



A step further go attem pts to outfit SPH with a Riemann solver (e.g. Ilnutsukal 120021 : 



Cha & WhitworthL 120031) . which makes SPH more similar to the Godunov approach 



used in Eulerian hydrodynamics. If successful, this may in principle eliminate the need 
for an artificial viscosity entirely, which would be a very welcome feature. So far, these 
versions of Godunov-based SPH have not found widespread application though, pre- 
sumably an indication that still more work is required to make them equally robust as 
the standard SPH formulation in all situations. 



6.3 Modelling Mixing 

As discussed earlier, the standard form of SPH has no mechanism to equilibrate the spe- 
cific entropies of neighboring particles, a property that does not allow different phases 
to mix into a perfectly homogeneous phase where all particle have the same thermody- 
namic properties. As seen above, this even causes an artificial surface tension effect. It 
has been suggested that the issue of mixing is ultimately the primary cause for the dif- 
ferences seen in the central entro py profiles of galaxy clusters simulated with SPH and 



mesh-based Eulerian techniques dMitchell et all 120091) . 



Recently, there hav e been a few suggestions to treat mixing in SPH by invoking arti- 
ficial heat conduction. Pried ( 20081) nas shown that this improves the results for Kelvin- 



Helmholtz instability tests when they are started with s harp discontinuities in the initial 
conditions. Also, in a study by lRosswog & Price d2007l) . it was found that this improves 



the accuracy of very strong shocks in Sedov-Taylor blast waves. One problem is to 
parameterize the heat conductivity such that it only affects the calcu l ations where it is 
really needed, and not anywhere else. The parameterization by IPricd (|2008|) . for exam- 
ple, also operates for static contact discontinuities, where it is not clear that any mixing 
should occur in such a case. A more physical parametrization of the diffusive mixing 
has been suggested by lWadsley. Veeravalli & Couchmanl ( 2008|) who use an estimate of 
the local velocity shear to parametrize the mixing rate. 

Yet another appro ach to mixing has been recently proposed bv Read. Hayfield & 
Agertz d2010h . who do not attempt to make the thermodynamic properties of a mixed 
fluid equal at the particle level. Instead, the number of smoothing neighbors is drastically 
enlarged by an order of magnitude, such that resolution elements of the mixed phases 
are sampled by a very large number of particles. 

6.4 Alternative Fluid Particle Models 

One possibility to substantially modify SPH lies in replacing the kernel interpolation 
technique of SPH by a different tec hnique to estima t e den sities. This idea has been ex- 
plored in the fluid particle model of iHeB & Springe where an auxiliary Voronoi 



tessellation was used to define a density estimate based on the mass of a particle divided 
by its associated Voronoi volume. The latter is simply the region of space that lies closer 
to the point than to any other particle. The same Lagrangian as in Eqn. (fT2l can then be 
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used to derive the equations of motion of the new particle based scheme, where now no 
smoothing kernel is necessary. This highlights that the kernel interpolation technique in 
principle only really enters in the density estimate of SPH. If the density estimate is re- 
placed with another construction, different particle-based hydrodynamical schemes can 
be constructed that are n evertheless similar in spirit to SPH. The Voronoi based scheme 
of iHeB & Spri ngel does not show the surface tension effect and has the advantage 

of leading to a consistent partitioning of the volume. It may also be more amenable to 
Riemann-solver based approaches to avoid the need of artificial viscosity, because an 
unambiguous interaction area between two particles (namely the area of the facet shared 
by the two Voronoi cells) is given. However, it is not clear whether it can significantly 
reduce the noise inherent in multidimensional SPH calculations. 

A considerable step further is the moving mesh technique devised by lSpringel (2009). 
This also exploits a Voronoi mesh, but solves the fluid equations with the finite-volume 
approach of the traditional Eulerian Godunov approach, where a Riemann solver is used 
to work out fluxes across boundaries. This is really a Eulerian scheme, but adjusted to 
work with particle-based fluid cells. In fact, as ISpringell (120091) discuss es, this approach 
is identical the well known Eulerian MUSCL scheme (Ivan Leer , 1984 ) when the mesh- 
generating points are kept fixed. If they are allowed to move with the flow, the scheme 
behaves in a Lagrangian way similarly to SPH. Nevertheless, in this mode the conceptual 
differences with respect to SPH also become most apparent. The first difference is that 
the errors in the discrete kernel sum are avoided through the use of an explicit second- 
order accurate spatial reconstruction. The second difference is that the cells/particles 
may exchange not only momentum, but also mass and specific entropy, allowing already 
for the mixing that may happen in multi-dimensional flow. 



7 CONCLUSIONS 



Smoothed particle hydrodynamics is a remarkably versatile approach to model gas dy- 
namics in astrophysical simulations. The ease with which it can provide a large dynamic 
range in spatial resolution and density, as well as an automatically adaptive resolution, 
are unmatched in Eulerian methods. At the same time, SPH has excellent conservation 
properties, not only for energy and linear momentum, but also for angular momentum. 
The latter is not automatically guaranteed in Eulerian codes, even though it is usually ful- 
filled at an acceptable level for well-resolved flows. When coupled to self-gravity, SPH 
conserves the total energy exactly, which is again not manifestly true in mesh-based 
approaches to hydrodynamics. Finally, SPH is Galilean-invariant and free of any er- 
rors from advection alone, which is another significant advantage compared to Eulerian 
approaches. 

Thanks to its completely mesh-free nature, SPH can easily deal with complicated 
geometric settings and large regions of space that are completely devoid of particles. 
Implementations of SPH in a numerical code tend to be comparatively simple and trans- 
parent. For example, it is readily possible to include passively advected scalars in SPH 
(for example chemical composition) in a straightforward and simple way. At the same 
time, the scheme is characterized by remarkable robustness. For example, negative den- 
sities or negative temperatures, sometimes a problem in mesh-based codes, can not occur 
in SPH by construction. Although shock waves are broadened in SPH, the properties of 
the post-shock flow are correct. Also, contact discontinuities can even be narrower than 
in mesh-based codes. 

All of these features make SPH a very interesting method for many astrophysical 
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problems. Indeed, a substantial fraction of the simulation progress made in the past two 
decades on understanding galaxy formation, star formation and planet formation has 
become possible thanks to SPH. 

The main disadvantage of SPH is clearly its limited accuracy in multi-dimensional 
flows. One source of noise originates in the approximation of local kernel interpolants 
through discrete sums over a small set of nearest neighbours. While in ID the conse- 
quences of this noise tend to be quite benign, particle motion in multiple dimensions 
has a much higher degree of freedom. Here the mutually repulsive forces of pressurized 
neighbouring particle pairs do not easily cancel in all dimensions simultaneously, espe- 
cially not given the errors of the discretized kernel interpolants. As a result, some 'jitter' 
in the particle motions readily develops, giving rise to velocity noise up to a few percent 
of the local sound speed. This noise in SPH is likely also a primary cause for the mixed 
results obtained thus far for MHD techniques implemented on top of SPH, despite the 
considerable effort invested to make them work accurately. 

In some sense it may seem surprising that SPH still models acceptable fluid behaviour 
despite the presence of this comparatively large noise. However, because SPH accurately 
respects the conservation laws of fluid dynamics as described by the Lagrangian, the 
local fluctuations tend to average out, enforcing the correct large-scale motion of the 
fluid. Fulfilling the conservation laws is hence clearly more important than a high order 
of the underlying scheme. However, as a consequence of this noise, individual SPH 
particles cannot be readily interpreted as precise tracers of the local state of the fluid 
at their representative location. Instead one needs to carry out some sort of averaging 
or binning procedure first. In this sense, SPH then indeed exhibits a kind of 'Monte- 
Carlo' character, as manifested also in its comparatively slow convergence rate for multi- 
dimensional flow. 

The noise inherent in SPH can be reduced by using a larger artificial viscosity and/or a 
higher number of smoothing neighbors. Both approaches do work at some level, but they 
are not without caveats and hence need to be used with caution. The standard artificial 
viscosity makes the simulated gas slightly viscous, which can affect the calculated solu- 
tions. The results will then in fact not converge to the solution expected for an inviscid 
gas, but to those of a slightly viscous gas, which is really described by the appropriate 
Navier-Stokes equations and not the Euler equations. However, improved schemes for 
artificial viscosity, such as a time-dependent artificial viscosity with judiciously chosen 
viscosity triggers, can improve on this substantially. Simply using a larger number of 
smoothing neighbors is computationally more costly and invokes the danger of suffering 
from artificial particle clumping. Here modified kernel shapes and different viscosity 
prescriptions may provide a satisfactory solution. 

In the future, it will be of primary importance to make further progress in understand- 
ing and improving the accuracy properties of SPH in order to stay competitive with the 
recently matured adaptive-mesh-refining and moving-mesh methodologies. If this can 
be achieved, the SPH technique is bound to remain one of the primary workhorses in 
computational astrophysics. 
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